version 7
#delimit;
set more off;
log using boehmke2006pa-graph.log, text replace;

/*	************************************************************************	*/
/*     	File Name:	boehmke2006pa-graph.do						*/
/*     	Date:   	May 12, 2006							*/
/*      Author: 	Frederick J. Boehmke						*/
/*      Purpose:	This file generates Figure 1, which shows the effect of		*/
/*			unobserved and observed characteristics on the probability	*/
/*			of voting for NAFTA with three separate scenarios.		*/
/* 			Observed factors are net endorsements. Unobserved factors are 	*/
/* 			accounted for by letting the observed position timing (timing)	*/
/* 			deviate from the expected timing at the mean values (timeavgD  	*/
/* 			and timeavgR).							*/
/*      Input File:	boehmke2006pa.dta						*/
/*      Output File:	boehmke2006pa-graph.log, boehmke2006pa-graph.dta		*/
/*			boehmke2006pa.gph						*/
/*	Version:	Stata 7.							*/
/*	************************************************************************	*/


	/********************************************************/
	/* This is the Weibull SUDCD likelihood.		*/
	/********************************************************/


program define expwbl;

	version 7;
	args lnf theta1 theta2 theta3 theta4;

	tempvar lambda1 lambda2 rhostar durdep P0_cond P_marg P_cum P_cens0 P_cens1;

	quietly gen double `lambda1'	= exp(-`theta1');
	quietly gen double `lambda2'	= exp(-`theta2');
	quietly gen double `rhostar' 	= (exp(2*`theta3')-1)/(exp(2*`theta3')+1);
	quietly gen double `durdep' 	= exp(`theta4');

	quietly gen double `P0_cond'	= 1 - exp(-`lambda1')*(1+`rhostar'*(2*exp(-($ML_y2*`lambda2')^`durdep')-1)*(exp(-`lambda1')-1));
	quietly gen double `P_marg'	= `lambda2'*`durdep'*((`lambda2'*$ML_y2)^(`durdep'-1))*exp(-(`lambda2'*$ML_y2)^`durdep');
	quietly gen double `P_cum'	= (1-exp(-`lambda1'))*(1-exp(-($ML_y2*`lambda2')^`durdep'))*(1+`rhostar'*exp(-($ML_y2*`lambda2')^`durdep'-`lambda1'));
	quietly gen double `P_cens0'	= (1-exp(-`lambda1')) - `P_cum';
	quietly gen double `P_cens1'	= 1 - (1-exp(-`lambda1')) - (1-exp(-($ML_y2*`lambda2')^`durdep')) + `P_cum';

	quietly replace `lnf' = ln(`P0_cond')   + ln(`P_marg') 	if $ML_y1==0 & rtcensr==0;
	quietly replace `lnf' = ln(1-`P0_cond') + ln(`P_marg') 	if $ML_y1==1 & rtcensr==0;
	quietly replace `lnf' = ln(`P_cens0') 			if $ML_y1==0 & rtcensr==1;
	quietly replace `lnf' = ln(`P_cens1') 			if $ML_y1==1 & rtcensr==1;

end;


	/********************************************************/
	/* Run the Weibull SUDCD model with parameterized	*/
	/* correlation. Code the right-censoring variable for	*/
	/* the likelihood function (for details, see 		*/
	/* comment in boehmke2006pa-analyze.do).		*/
	/********************************************************/

use boehmke2006pa;

  generat _rtcens = rtcensr;

ml model lf expwbl (vote= contdiff mexbordr pscenter hhcenter partyid numdiff) 
	(timing=corptpct labtpct mexbordr dleader rleader ncomact ideol pscenter inter1 hhcenter inter2) 
	(rho: partyid) /durdep;

	ml search;
	  ml init rho:_cons=0;
	  ml init rho:partyid=0;
	ml maximize, difficult;

	/********************************************************/
	/* Set up the data to generate the predicted values. 	*/
	/* After getting the means, expand to create 		*/
	/* observations that vary the timing of announcement.	*/
	/* There are three sets of predicted values, each 	*/
	/* calculated for 500 days and indexed by simno.	*/
	/* 							*/
	/* 1. Dems - varies unobserved factors only.		*/
	/* 2. Dems - varies observed and unobserved factors.	*/
	/* 3. Reps - varies observed factors only.		*/
	/********************************************************/


		/********************************************************/
		/* First for Republican legislators (partyid=0).	*/
		/********************************************************/

collapse (mean) corptpct labtpct pscenter hhcenter numdiff contdiff 
	(p50) mexbordr dleader rleader ncomact ideol partyid if partyid==0;

  expand 500;

  generat timing = _n;
  generat simno = 3;

  save boehmke2006pa-graph, replace;


		/********************************************************/
		/* Then for Democratic legislators (partyid=1).	They	*/
		/* have two version - with observed endorsements 	*/
		/* (simno=2) and without (simno=1).			*/
		/********************************************************/

use boehmke2006pa;

  collapse (mean) corptpct labtpct pscenter hhcenter numdiff contdiff
	(p50) mexbordr dleader rleader ncomact ideol partyid if partyid==1;

  expand 1000;

  generat simno = 1 if _n <= 500;
  replace simno = 2 if _n >  500;

  generat timing = _n;
  replace timing = _n - 500 if _n > 500;

  append using boehmke2006pa-graph;
  save boehmke2006pa-graph, replace;


	/********************************************************/
	/* Now add in the observed position announcement data	*/
	/* to the just generated data of mean values. Use 	*/
	/* replace command to fill in values for days not in	*/
	/* the original data set. Then set them to zero for the */
	/* correlation-only simulation for Democrats (simno=1). */
	/********************************************************/

use boehmke2006pa;

  keep numdiff timing;

	collapse numdiff, by(timing);
	joinby timing using boehmke2006pa-graph, unmatched(using);

  sort simno timing;

  	replace numdiff  = numdiff[_n-1]  if numdiff  == . & _n != 1 & simno==simno[_n-1];
  	replace numdiff  = 0 if simno==1;

  save boehmke2006pa-graph, replace;

	/********************************************************/
	/* Recode interactions to product of mean values.	*/
	/********************************************************/

  generat inter1=ideol*pscenter;
  generat inter2=ideol*hhcenter;

	/********************************************************/
	/* Generate predicted values of the indices (xb) and	*/
	/* transform them into intermediate variables to use 	*/
	/* to calculate predicted durations at mean values. 	*/
	/********************************************************/

  predict votepred_xb, eq(#1) xb;
  generat lambda1 = exp(-votepred_xb);

  predict timepred_xb, eq(#2) xb;
  generat lambda2 = exp(-timepred_xb);

  generat rhostar = (exp(2*[rho]_b[_cons] + [rho]_b[partyid]*partyid) - 1)/
		    (exp(2*[rho]_b[_cons] + [rho]_b[partyid]*partyid) + 1);

  local p = exp([durdep]_b[_cons]);


	/********************************************************/
	/* Now generate the predicted durations and votes at  	*/
	/* the mean values. First generate local scalar 	*/
	/* variables to insert lines in the graph command. 	*/
	/********************************************************/

local timeavgD 	= round(exp(lngamma(1 + (1/`p')))*exp(timepred[1]),1);
local timeavgR 	= round(exp(lngamma(1 + (1/`p')))*exp(timepred[1001]),1);

  display "Average Democratic Announcement: " `timeavgD';
  display "Average Republican Announcement: " `timeavgR';

local voteavgD = exp(-lambda1[1])*(1+rhostar[1]*(2*exp(-(`timeavgD'*lambda2[1])^`p')-1)
	*(exp(-lambda1[1])-1));
local voteavgR = exp(-lambda1[1001])*(1+rhostar[1001]*(2*exp(-(`timeavgR'*lambda2[1001])^`p')-1)
	*(exp(-lambda1[1001])-1));

  display "Average Democratic Vote Probability: " `voteavgD';
  display "Average Republican Vote Probability: " `voteavgR';


	/********************************************************/
	/* Finally, generate the predicted vote probabilities	*/
	/* given that the actual announcement time deviates 	*/
	/* from the expected time at the mean values. Put it 	*/
	/* into three separate variables for graphing.		*/
	/********************************************************/

generat votepred = exp(-lambda1)*(1+rhostar*(2*exp(-(timing*lambda2)^`p')-1)*(exp(-lambda1)-1)) 
	if simno==1 | simno==2;
replace votepred = exp(-lambda1) if simno==3;

generat voteprD1	= votepred if simno==1;
generat voteprD2	= votepred if simno==2;
generat voteprR		= votepred if simno==3;


  label variable timing 	"Timing of Position Announcement";
  label variable votepred 	"Probability of Pro-NAFTA Vote";

  label variable voteprR  	"Republicans (Observed Only)";
  label variable voteprD1 	"Dems (Unobserved)";
  label variable voteprD2 	"Dems (Unobserved and Observed)";

  save boehmke2006pa-graph, replace;


	/********************************************************/
	/* Now generate the graph. Only use every fifth day 	*/
	/* (via the if command) to make the graph more readable.*/
	/********************************************************/

twoway connected voteprR voteprD1 voteprD2 timing if int(timing/5)==timing/5, scheme(s1mono)

	msymbol(O i Th) msize(1.3 1.3 1.3) lpattern(blank solid blank) lwidth(medium medthick medium)
	mcolor(gs0 gs0 gs0)

	xline(`timeavgD') yline(`voteavgD' `voteavgR') 
	ylabel(0.25 0.5 0.75) xlabel(0 50 100 150 200 250 300 350 `timeavgD' 463 500)
	ytitle("Probability of pro-NAFTA position")

	text(0.55 `timeavgD' "Expected Announcement Day", place(w)) 
	text(`voteavgD' 0 "Baseline Democratic Probability", place(ne)) 
	text(`voteavgR' 0 "Baseline Republican Probability", place(se)) 

	legend(cols(1) rows(3) order(2 3 1))

	saving(boehmke2006pa, replace);

log close;
clear;
exit;
